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Multigrid Methods for Space Fractional Partial Differential 

Equations* 

Yingjun Jiang^and Xuejun Xu^ 


Abstract 

We propose some multigrid methods for solving the algebraic systems resulting from finite 
element approximations of space fractional partial differential equations (SFPDEs). It is shown 
that our multigrid methods are optimal, which means the convergence rates of the methods are 
independent of the mesh size and mesh level. Moreover, our theoretical analysis and convergence 
results do not require regularity assumptions of the model problems. Numerical results are given 
to support our theoretical findings. 
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1 Introduction 


Fractional partial differential equations (FPDEs) have found many impressive applications in lots 
of fields, such as finance, phase transitions, stratified materials, anomalous diffusions (see [25] and 
references therein). To solve them, both analytical and numerical methods are used in the literature. 
The analytical methods like the Fourier transform method, the Laplace transform method and 
the Mellin transform method have been developed to seek closed-form analytical solutions |32j . 
Since such closed-form analytical solutions are unavailable in most cases, extensive researches have 
already been carried out on the development of numerical methods for fractional partial differential 
equations like finite difference methods (see e.g., |5l[TTl[T8l[26l[271|36ll39]), finite element methods 
(see e.g., [T2lfT^[23] L and spectral methods [191121] . 

Let 12 be a polyhedral domain in we consider the space fractional partial differential equa¬ 
tions (SFPDEs): find u{x) such that (see [T5] i 


Isd- 


D‘1°‘u{x) M{z)dz + cu{x) = f{x) 


X 


G D, 


( 1 . 1 ) 


'^\R‘i\Q — 0 ) ( 1 - 2 ) 

where 1/2 < a < 1, c> 0, / is a source term, || • ||2 denotes the standard Euclidean norm, 
= {z € \\z \\2 = 1}, M{z) is a probability density function on and which will 

be given later, denotes the directional derivative of order 2a in the direction of the unit vector z. 
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Here we assume M is symmetric about origin, i.e., M{z) = M{z') if z, z' G satisfy z + z' = 0, 
which means that the considered problem is a symmetric one. 

One special case of a is 


d 

-^{Pi -ooDI"^ + qi XiD‘^^)u + CU = f (1.3) 

i=l 

and Pi,qi > 0 satisfying pi = qi and Yli=i{Pi + 9*) = where -ooD'i^, xiD'^ denote Riemann- 
Liouville fractional derivatives. Actually, (11.31) can be obtained from (II.ip by taking M = Yli=i Pi^{^~ 
Si) + qiS{z + Si), where Cj is the ith column of identity matrix in and 6 the Dirac function 

on The corresponding time-dependent equation of (II.ip can be used to describe a general 

super-diffusion process (see [21])) which is an appropriate extension from one dimensional problem 

I 

— - (p -ooDI^ + q xD^)u + cu = f. (1.4) 

As to the super-diffusion, please refer to |2S] for details. 

One of the greatest challenges for numerically solving SFPDEs is how to reduce the computation 
costs. Due to the nonlocal properties of fractional differential operators, numerical methods for 
linear SFPDEs tend to yield the linear equations Ax = b with the following characteristics: 1). 
the coefficient matrix A is dense or full; 2). the condition number of A increases fast, as the 
mesh becomes fine. Reducing the computation costs for SFPDEs is harder than doing it for the 
integer order PDFs. Some methods have already been designed to overcome this difficulty, such as 
alternating-direction implicit methods (ADI) |271ll2lll3|, and iterative methods [20113111431146] . 

Iterative methods seem to be efficient tools for solving SFPDEs. Actually two issues in this sit¬ 
uation need to be concerned for efficiency: one is to do the matrix-vector multiplications efficiently, 
and the other is to find good preconditioners. As to the first issue, some literatures are contributed: 
in with the notice of Toeplitz-like structure of the coefficient matrix, the matrix-vector mul¬ 
tiplications are done with 0(iV log A^) complexity by using a fast Fourier transform (FFT) |8l|9]. 
This technique of ” matrix-vector multiplication” has been widely used to improve the efficiency 
of iterative methods for the SFPDEs [2ni[3n 1431146] . As regards the second issue, some literatures 
should be listed as follows: the first relevant paper may be [2] in which a multilevel preconditioner 
of fractional power was put forward; in [20|, the authors propose preconditioners constructed by 
some banded matrices of fixed band width; in [45], the authors present a preconditioner by some 
symmetric positive Toeplitz matrixs; moreover a new preconditioner is designed in m through 
some circulant matrixs. 

It is known that multigrid methods are optimal iterative procedures, which have been widely 
used for integer order PDFs (see e.g., [3ll38] l. In recent years, some researchers begin to investigate 
multigrid methods for solving SFPDEs. For instance, in |4S], Zhou and Wu apply the multigrid 
method to solve one dimensional steady SFPDEs, and in [31] . the authors consider the V-cycle 
multigrid method for solving corresponding time-dependent problems. But till now, no satisfactory 
convergence results have been obtained for the multigrid methods for solving SFPDEs. Actually, 
in m, the authors only conduct the theoretical analysis for the two-level multi-grid method, and 
Zhou and Wu in |49| get the convergence results only under the assumption that the adjoint problem 
hold sufficiently smooth solution. 

In this paper, we introduce a V-cycle multigrid method with one smoothing step on each level 
to solve linear algebraic systems resulting from the finite element approximations of the SFPDEs 
dni). It is shown that our V-cycle multigrid methods are optimal, which means the convergence 
rates are independent of the mesh size and mesh level. Moreover, our theoretical analysis and the 


2 


convergence results in this paper do not require any regularity assumptions of the model problems. 
To the best of our knowledge, this paper is a first attempt to give a rigorous theoretical analysis for 
the V-cycle multigrid methods for the finite element approximations of SFPDEs in any dimensions. 

This paper is also the first work to design the fast solver for the SFPDE (|1.1|) with M being 
a continuous function. Among the current numerical methods for SFPDEs, most of them are for 
one dimensional problems and for some special high dimensional problems like (II.3p . and only a 
few are for more general problems like m- Actually, only [T6l[33] study the numerical methods 
for (|l.ip : in [16], the authors consider the finite element approximation for (|l.ip and in [33|, the 
author studies the corresponding time-dependent case. 

In the rest of the paper, no loss of generality, we restrict ourselves to the case d = 2, namely, we 
consider the problem in M^. Eor A C denote T^(A) the space of all measurable function v 
on A satisfying J^{v{x))‘^dx < oo, and C^PA) the space of infinitely differentiable functions with 
compact support in A. Set 

{v,w)a= / vwdxdy, ||u||a = 

Ja 

and they are abbreviated as {v,w) and ||u|| respectively if A = R^. 

To simplify our statement, we make a convention here: function v defined on a domain A C R^ 
also denotes its extension on R^ which extends v by zero outside A. The constant C with or 
without subscript will denote a generic positive constant which may take on different values in 
different places. These constants will always be independent of the mesh sizes and levels in the 
multigrid methods. Following m. we also use symbols <,> and ss in this paper. That ai < 6 i, 
“2 ^ ^2 and 03 PS 63 mean that oi < Ci^i, 02 > 6*262 and 6*363 < 03 < 6*363 for some positives 
61 , 62,03 and 6 '. 

The rest of the paper is organized as follows: for the sake of completeness, in section ITTI we give 
our model problem and the corresponding finite element discretization. In section [S] we present 
our V-cycle multigrid methods and introduce some basic theoretical results. In section 01 we shall 
prove the convergence of the multigrid methods. Finally in section 5, the numerical results are 
given to verify our theoretical findings. 


2 The model problem and its discretization 

In this section, we shall present the SFPDE in R^, and then introduce its variational formulation 
and corresponding finite element discretization. 


2.1 The model problem 


We first introduce the concepts of directional integrals and derivatives [TB] . 

Definition 2.1. [T 6 | Let /r > 0, 0 € R. The /ith order fractional integral in the direction z = 
(cos 0 , sin 0 ) is defined by 

/■oo ^//-l 

D~^v{x,y) := Dg'^v{x,y) = / v{x — t cos6,y — t sm6)dT, 

Jo 


where T is the Gamma function. 


Definition 2.2. [T 6 | Let n be a positive integer, and 0 G R. The nth order derivative in the 
direction of z = (cos 0 , sin 0 ) is given by 


Dev{x,y) 


„ d . a ^ 
cos —h sm y— 
ox oy 


n 


v{x,y). 
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Definition 2.3. [16] Let /i > 0, 0 E M. Let n be the integer such that n — 1 < fj, < n, and define 
a = n — fj,. Then the /rth order directional derivative in the direction of z = (cos 6, sin0) is defined 
by 

Di^v{x,y) := = D^D-^v{x,y). 

If V is viewed as a function in x, Dq, are just the left and the right Riemamm-Liouville 
derivatives (see e.g., [32l[35]). The fractional derivative operators in problem (ini) are related to 
the following fractional derivative: 

Definition 2.4. |16] Assume that v : M? ^ W, y > 0. The yth order fractional derivative with 
respect to the measure M is defined as 


D'"^'>^{x,y) := f D^v{x,y)M{e)de, 

where = [0 + v, ‘I'k + v) with a suitable scalar and M{9), which satisfies M{9)d9 = 1, is 

a periodic function with period 27r. Usually we take = 0, if it causes no unreasonable expression 
(see (|2.2p L 


Remark 2.5. It is easy to check that 


Dljv{x,y) = an 


d'^v d'^v d'^v 

TTW + fl22 7rw + 2ai2W— 
ox^ oy^ oxoy 


where an = cos"^ 9M(9)d9, 022 = s\v? 9M{9)d9 and ai 2 = 2 cos 9 sin 9M(9)d9 (see 
also [25]). Denote L a positive integer, let 9^ E [0,27r) and pk > 0, k = 1,2,...,L, satisfy 
Ylik=iPk = 1- Assume that D^v is continuous in 9, and then 

L 

= ( 2 . 1 ) 

k=l 

if 

L 

M = Y,Pk5{9-9u), ( 2 . 2 ) 

k=l 

where 5 denotes Dirac delta function. 

For tt : ^ M, define differential operator in as 

LaU = —DI^u + cu. 


Denote Q a polygonal domain in set 1/2 < a < 1, and then the model problem of this paper 
is to find u : D —>■ M such that 

f LaU = f, in n, , . 

\ a = 0, ondn, ^ ’ 

where / is a source term and we assume that M{9) satisfies M{9) = M{9 + vr) for 0 E M, i.e., (|2.3I) 
is a symmetric problem. Here, we recall the convection made in Section (TJ i.e., u also denotes its 
extension by zero outside D. 
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2.2 The variational formulation 


Definition 2.6. [30] Let /U > 0, be the Fourier transform of v{x,y), |^| = \/ + ^ 2 - 

Define norm 

(1 + I^IWH • 

Let := {u G L^(]R^); Hul< oo}. 

For V G ifg(D), we also denote ||u||j:^m(r 2 ) by ||u||//M(n)- It is known that is a Hilbert 

space equipped with the inner product (u, r(;)j:/M(R 2 ) = ((1 + \C\'^)^Tv,Tw) and is dense in 

(see m)- Now, we introduce and prove some useful results for the fractional directional 
derivatives of functions in 

Lemma 2.7. For /r G M, u G the Fourier transform of DqV is 


\v\\hi^ 


F{D'^v{x,y)) = (27rz(^icos(9 + ^2sin(9))'" 
Lemma 2.8. For ii,s > 0, v,w £ 


{D^v,w) = {D^ %,Dl^^w), 

where D^v = v. 

Proof. By lemma ITTl and (|A.1|1 . we know FDgV = (27rz(^i cos0 + ^2 sin0))^ FDq~^v = 
(27ri(,fi cos 0 + ^2 sin0))^~* FDqj^^w = (27ri(^i cos0 + .^2 sin0))®J^t(;. Then the lemma follows 
by Parseval’s formula. □ 

We define the weak fractional directional derivative according to the relation {DqV,w) = 
{v,Dgj_^w) which is a special case of Lemma 12.81 (see also Lemma 5.7 in [l6|). Let 
denote the set of locally integrable functions on 

Definition 2.9. Given // > 0, 0 G M, let u G L^(M^). If there is a function G such that 

(^^, Dg+.^w) = {v^,w), Vu; G ^^(K^), 

then is called the weak yih. order derivative in the direction of 0 for u, denoted by D^v, i.e., 
Vf, = D^v. 

It is not hard to see that the weak derivative DgV is unique if it exists and that the weak 
derivative coincides with the correspondent derivative defined in Definition 12..SI if v G C'^(M^). In 
the following, we use D^v to denote the weak derivative. 

Lemma 2.10. Let y > 0. For any v G 0 < s < y and 0 G M, the weak derivative DgV 

exists and satisfies 


J^Dgv{^i,^2) = (27rf^icos0 + 27rz^2sin0)^T'u(6,6), (2-4) 

ll-D^ull < C'||u||j^^(]r2). (2.5) 

Proof. Since C'“(M^) is dense in F7^(M^), there is a Cauchy sequence {vn} C C'“(M^) such that 
Ikn — '*^||hm(r 2 ) —>■ 0 as n —0. By lemma ITTl FDqW = (27rz(^i cos 0 + ^2 sin 0))* T'rc for w G 
C'“(M^). By Parseval’s formula and 0 < s < /i, it is not hard to see that ||Dgt(;|| = ||J-'D|tc|| < 
C'||t(;||j:^M(R 2 ). So we have ||D|u„ — DgVmW ^ C\\vn — Vm\\Hr{R^) and {DgVn} is a Cauchy sequence 
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in Denote Vs E the function to which {D^Vn} converges to. By Lemma 12.81 for any 

w G Co“(M2 ), 

{Vn,Dg^^w) = {DgVn,w). 

Take the limits of both sides of the above equation, we obtain {v, Dg^^w) = {vs,w) for any w G 
So DgV exists and is equal to Vg by Definition 12.91 By the definition of Fourier transform 
for the function in 

{i2Tri{^icos9 + ^2sme)YTvn,v) = iDlvn,J^v), VuGCo“(m2). (2.6) 


Because 


v„. - V 


Hi^ 


{l + \erl-^\HVn-v)\ 


0 , 


it is not hard to see that (27ri(^i cos 0 + ^2 sin 9)Y Fvn converges to (27ri(|i cos 0 + ^2 sin 9)Y Fv in 
Take the limits of both sides of (|2.6H . we obtain (|2.4I) by the definition of Fourier transform. 
(12.51) can be proved directly by (|2.4I) and Parseval’s formula. □ 


Lemma 2.11. Let fj,,s > 0 with fj, — s > 0. For v,w G 


iD^v,D^^^w) = {D^+^v,D^Y> 


(2.7) 


Proof. For any g G 7 7^+^(R ^), \\D^g\\, \\D^j^^g\\, g\\ and are all bounded by 

by Lemma I2.1UI Then the lemma follows from that C'^(R^) is dense in and 

Lemma 12.81 □ 

Assume that the solution u of (j2.3p is sufficiently smooth (indeed, that u G C'^(D) with u|gn = 0 
is sufficient). Multiplying both sides of the first equation in (12.31) with v G (7^(11) and integrating 
over D give 

{Dg°‘u,v)M{9)d9 + c{u,v) = {f,v), v € C^(fl). (2.8) 

Then employing the relation {DgW,v) = {w,Dl_^^v) (it can be obtained by integration by parts), 
we obtain 

{Dl^-\Dl^,v)M{9)d9 + ciu,v) = (/,u), v G C^{n). (2.9) 

Then by Lemma 12.111 (|2.9p can be rewritten as 

(D>, D^^^v)M{9)d9 + c{u, v) = (/, u), vG C^{n). (2.10) 





Define the bilinear form B : L7g (D) x L7 q (D) —)• R as 

B{u,v) := - / {DgU,Dg_^_.j^v)M{9)d9 + c{u,v). 

Jo 

By M{9) = M{9 + vr) for 0 G R, it is easy to check that B{v,w) is a symmetric bilinear form, i.e., 
B{v,w) = B{w,v) for v,w G Hq{Q). The variational formulation of (12.31) is (see also [16]) to find 
u G Hq{LI) such that 

B{u,v) = if,v), VuGFo“(D). (2.11) 

Now we restate some results in m about the solvability of (|2.11l) . To guarantee the existence of 
the solution of (12.lip , we assume that M{9) satisfies 

p27r 

/ |(eicos 0 + 6 sin 0 )| 2 “M( 0 )d 0 >Co|eP“ ( 2 . 12 ) 

Jo 
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for some positive Cq. Denote k = 27r(^i cos0 + ^ 2 sin0), Ei = {(Cij^ 2 ) E : ^1 cos0 + ^ 2 sm0 > 0}, 
E 2 = {(Cl) ^ 2 ) € : Cl cos 0 + C 2 sin0 < 0}, and then by Parseval’s formula and Lemma [2.101 


{Dfv, 


m+.v) 


(|k|^“ exp(iasign(K)7r)J^i;, Fv) 

(|kP“ exp(ia7r)J^i;, Fv)ei + exp(—ia7r)J^i;, Fv)e 2 
cos(Q!7r)(|Kp“J^u, Fv) 

+isin(a7r) [{\k\^^Fv,Fv)e^ — {\h\^^Fv,Fv)e 2 ) 

cos(Q!7r)(|Kp"J^i;, J^u), (2.13) 


where the computation of complex please refer to Appendix, in the fourth equality, the Euler 
formula exp(m) = cos(k) + i sin(K) is used, the last equality is because the value of {DgV, Dg^^v) 
is real and the imaginary part must be zero (another proof for this equality please refer to |16ji. 
Furthermore, by (j2.12p and cos(a7r) < 0 


> 

rv-/ 


/*27r 

/ iD^v,D^^^v)M{e)de 

Jo 

r r /*27r 

(ott) // \Fv\^ / |27r(Ci COS0 + C2sin0)p"M(0)(i0(iCif^C2 

J JR2 Jo 

\i\^^\Fv\^di^di2. 


— cos 


(2.14) 


For V G Hq{Q), we have 

\\vf < Ci\\D^v\\^ = Ciff |27r(Cicos0 + C2sin0)|2“|.Fu|2dCidC2 

J Jr 2 

< CI 2 // |CP“|-Fi;pdCidC2, (2.15) 

J Jr 2 

where the inequality is by (5.15) in [16] and the equality is by Parseval’s formula. With the 
combination of ([2.1411 and ([2.1511 . we conclude under condition ([2.1211 . 

B{v,v)>\\v\\jjc,^^^, veH^in). (2.16) 

By Lemma 12.101 it is easy to verify that 

B{v,w) v,weHo{n). (2.17) 

By (I2.16P and (I2.17p . using Lax-Milgram theorem, we know that the variational formulation ([2.1111 
admits a unique solution in 77g (D). 

Remark 2.12. Condition (I2.12p is easily satisfied. For example, it holds if M{9) is non-zero over 
a connected set of positive measure in [0,27r) (see [E]), and it holds when M{6) = Yl^=iPk^{^ ~ 
k'K/2)d6, with pk > 0 and Pi + Ps = 1, P 2 F P4 = 1- 


2.3 The finite element discretization 

Let 7)i be a quasi-uniform triangulation of D such that D = Ui^-gT^iP, hx be the maximal length 
of the sides of the triangle K and h = max^e?^ hx- Denote Pi{K), I > 1, the space of polynomials 
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of degree less than or equal to I on K £ Th- Define the finite dimensional subspace V associated 
with Th as 

D := {u G C°(D) : u|an = 0,v\k G Pi{K)yK G Th}. 

It is known that V C Hq{Q,) C Hq{Q). Thus the finite element approximation for ( 12 . 111 ) is to find 
Uh £V such that 

B{uh,v) = if,v), VuGD. (2.18) 

The error estimates for the finite element solution Uh are given in m- 

In practical applications, we use the finite element discretization (j2.18ll only when the probabil¬ 
ity density function M has the discrete form as that in (12.2p (when M {9) is a continuous function, 
the finite element discretization (I2.18P can hardly be realized). For the case that M{9) is the con¬ 
tinuous function, we propose an alternative finite element discretization instead of (j2.18p . Here we 
focus on the case M{9) G 6*^(0, 27r] is a periodic function with period 27r to present our alternative 
finite element problem: find Uh £V such that 

B{uh,v) = {f,v), VuGH, (2.19) 

where B{.,.) is an approximation of i?(-,-). Exactly in this paper, set a positive integer Nq such 
that Nq is a multiple of 4. Letting 9i = 2171 /Nq, z = 0,..., — 1 and denoting A0 = 2ttINq, we 

use the compound trapezoid formula to get B{., •), i.e., for v,w gV, 

/*27r 

B{v,w) = - {DQV,Dg_^_^w)M{9)d9 + c{v,w) 

Jo 

Ng-l 

- {Df.v, + c{v, w) := B{v, w). 

i=0 

The fact that M(9) = M(9 T tt) and Ng is a multiple of 4 guarantees that B{v,w) is a symmetric 
bilinear form as well, i.e., B{v,w) = B{w,v). By Parseval’s formula, we have 

{DgV,Dg^^w)n = {{2Tri^i cos 9 + 2Tri^2 sin 9)'^'^Tv, Tw) 

< C||u||//«(n)||u;||//a(f2) (2.20) 


and 

d 


^^{DgV, DQ_^_^w)n = 2q; ((— 27rz^i sin0-|-27rz,f2 cos0)(27rz^i cos 0-|-27rz^2 sin0)^“ ^Tv,Tw) 

< C||u||//a(f2)||u>||//a(f2). (2-21) 


By the error formula for the compound trapezoid formula, it is easy to verify that 

\B{v,w) - B{v,w)\ < C'A6»||u||H«(n)||w^||H“(f7), 


( 2 . 22 ) 


where C is a positive constant independent of 9, v and w. Combining (12.221) with (I2.16D and (I2.17D . 
we know for sufficiently small A9, 


B{v,w) < ||u||j^«(n)||u;||H-(n), v,w G LIq (^)- 


(2.23) 
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By Lax-Milgram theorem, (I2.19P has a unique solution. The first Strang lemma (see |10] i holds 
here, i.e.. 


< 


■ugv ^ ,, , 

C* { 11 ^^ - + A 6 »||u||j^«(n)} 


-II ^^-rlii II - B(v,w)\ \ 

U-Uh\\H<-m ^ C” inf S P - I' LH-“(n) + sup -- ) 

I w&v in b“(n) J 


Finally, the finite element approximation of (|2.3I) is unitedly presented as: find Uh G V such 
that 

B{uh,v) = if,v), VvGV, (2.24) 

where B{v, w) = — {DgV, DQj^^w)M{6)d6+c{v, w), M{9) is equal to a discrete form Yl,k=iPk^i^~ 
9k) such that B{-, •) is a symmetric bilinear form, 

B{v,v) > \\v\\ho,(^q),B{v,w) < ||u||/^«(q)||u;||h-(o), v,w G Fq (^)> (2-25) 

and M{9)d9 < 1. Specially for the cases mentioned above, the finite element problem (|2.24l) 

represents problem ()2.18p if M{9) = M{9) = Ylk=iPk^{^ — ^fc) and problem (|2.19p if M{9) = 
- 9i)M{9i). 


3 Multigrid algorithm 

In this section, for (|2.24ll . we shall present our V-cycle multigrid algorithm and a general framework 
for our convergence analysis. 

Take fuGV such that {fh,v) = (/, u), Mv G V and define a linear operator A : V ^ V bs 
follows: 

{Av,w) = B{v,w), 'iv,wGV. (3.1) 

The finite element approximation of system (I2.24p can be restated as to find Uh GV such that 

Auh = fh- (3.2) 

In the following, we shall use the operator equation (13.21) to construct our multigrid algorithm. 
Since B{v,w) is a symmetric bilinear form, we know, by (I2.25p . that A : V ^ V is symmetric 
positive definite with respect to (•, •), i.e., 

{Av,w) = {v,Aw), v,wGV-, {Av,v) > 0, O^vGV. 


Then bilinear form 

{v,w)a ■= {Av,w), v,wGV, 
also induces an inner product on V. Set norm 

\\v\\a = [Av^v)^^"^, vGV. 

By (|2.25l) . we have 

ll'^^IU ~ Vu G y. 


(3.3) 
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3.1 Algorithm 

Assume that the triangulation 7/i of 0 is constructed by a successive refinement process. To be 

precise, let 7j = 7h for some J > 1, and Tfc for A: > 0 be a nested sequence of quasi-uniform 

triangulations, i.e., Tk = {t^} consists of simplexes of size hk such that 12 = UjT^; is a union 
of simplexes of r^. We further assume that there is a positive constant 7 < 1, independent of k, 

such that hk is proportional to 7 ^ and the simplexes in 71 are of diameter ~ 1 . 

For each partition 71, we may define finite element spaces I 4 by 

Vfc = {u G C%n) : u|an = 0,u|, G P/(r),VT G Tk}. (3.4) 

Obviously, the following inclusion relation holds: lA C V 2 C • • • C Vj = 1/. Our V-cycle multigrid 

methods are based on the subspace decomposition V = lA -|- VI Vj. 

For each /c G {1, 2,... , J}, define projectors Qk, Pk ■ V ^ Vk by 

{QkV, w) = {v, w), {PkV, w)a = (v, w)a, V eV,w eVk, 

specially, set Qo : A as Qqv = 0, and define the linear operator Ak '■ Vk ^ 14 

{AkV,w) = {Av,w), v,weVk. 

It is easy to verify that 

AkPk = QkA, k = l,2,...,J. (3.5) 

It is obvious that Ak is symmetric and positive definite with respect to (•,•). Denote Xk G M, 
k = 1, 2,... , J, the maximal eigenvalue of Ak- 

Let Uk = Pk^h and fk = Qkfh: we may get the operator equation in subspace 

= fk- (3.6) 

Our multigrid algorithm is essentially an iterative procedure in which the subspace equation ()3.6I) 
is approximately solved successively to get new approximations to (|3.2p from old approximations. 
More precisely, denote : 14 —>■ I 4 the approximate inverse of Ak , and u°^‘^ the old approximation 
to u. Correcting the residual of in I 4 gives 

+ RkQk{fh - AT^^). 

We take Rk to be symmetric with respect to (•, •) such that 

{RkV,v) ^ Vu G 14,A: = 1,2,... , J. (3.7) 

Afc 

Remark 3.1. In this paper, we have hi = 0(1) and take i?i = By Lemma 14.31 ()3.3p and the 
definition of norm || • we know that {v,v) < {Aiv,v) < /ij“^“(u,u), and Ai = 0(1). Then 

we have (Aj“^u,u) ps ;^j^(u,u). 

Next we give our V-cycle multigrid algorithm. 

V-cycle Multigrid Algorithm. Let = 0 G V, assume that £ V has been obtained. 
Then is generated by 

u>^+^ = u^ + Bj{fk-Au^), (3.8) 

where Bj is defined inductively: Let Bi = and assume that Bk-i : 14_i Vk-i has been 
defined; then for g £ Vk, Bk '■ Vk ^ Vk is defined as follows: 

Step 1. = Rkg] 

Step 2. = v^ + Bk-iQk-iig - A^u^); 

Step 3. Bkg = v'^ + Rkig - AkV^). 
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3.2 A general framework 

For the V-cycle multigrid method, we have 

= {I-BjA){uh-u^). 


Denote 


Ej = {I-Tj){I-Tj_i)---{I-Ti), E*j = {I-T,)---{I-Tj_i){I-Tj) (3.9) 


with Ti = Pi, Tk 
norm as 


RkA^Pk, k = 2,3,..., J. Then we have {I — BjA) 


\\Ej\\a = sup 

v&V 


IblU 


It is easy to see that Ej is the (•, •)A-adjoint of Ej, i.e., 


EjEj. Define the operator 


{Ejv, w)a = {v, E*jw)a, v,w £V 


and that 

||F;j|U = ||i^}IU, \\EjE*j\\a<\\Ej\\\. 

The main work in this paper is to establish the contraction property: there is a constant 
0 < 5 < 1 independent of the mesh size and mesh level such that 

||^j|U<V^. (3.10) 


By (|3.10p . we may obtain \\uh — u^\\a < S^\\uh — 

Remark 3.2. For the V-cycle multigrid method, the spectral radius of the iterative matrix p = 
p{I — BjA) < S. It is known that the condition number k{BjA) < and BjA is self- 

adjoint and positive with respect to inner product (•,-)^. The d’s independence of the mesh size 
implies that Bj is a good preconditioner for A which can be used to design efficient preconditioned 
conjugate gradient methods. 

Define Kq and Ki as two smallest positive constants satisfying the following conditions: 

1. For any u G V, there exists a decomposition v = Xlii for Vi G Vi such that 

J 

^{R~^Vi,Vi) < Ko{Av,v). (3.11) 

i=l 

2. For any S' C {1, 2,... , J} x {1, 2,... , J} and Vi,Wi G V for i = 1, 2,..., J, 



(J Y 

/ J 


{TiVi,TjWj)A < Ki 

'^{TiVi,Vi)A 

\'^{TjWj,Wj)A\ ■ 

(3.12) 

(*d)e5 

Vi=l / 

\3=l ) 



The estimate of the upper bound of HF'jm relies on the following lemma: 
Lemma 3.3. Let Ej he defined by We have 




2 — uji 

Ko{l + Ki)2' 


where oji = va&y.p{RkAk), p{RkAk) denotes the spectral radius of RkAk- 
k 
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The estimate of the parameter coi is straightforward. 
(|3.7p . foT V eVk {k = 2,..., J) 


Since Ri = p{RiAi) = 1 . 


From 


^iv,v) < {RkV,v) < 


and furthermore 


C 2 

{RkAkV,v)A = {RkAkV,Akv) < —{AkV,Akv) < C2{v,Akv) = {v,v)a-, 

Afc 


(3.13) 


where the last inequality is obtained from that Ak is symmetric positive matrix and \k is the 
maximal eigenvalue of Ak- Combining (|3.13p with the fact that RkAk is symmetric with respect 
to inner product {■,-)a, we have p{RkAk) < 6 * 2 . Taking Rk such that C 2 is suitably small can 
guarantee the wi < 2 . 

Next, we shall estimate the parameters iFi, K 2 . The following Lemma is helpful for the analysis. 

Lemma 3.4. Let e = {eij) G R^^^ he a nonnegative symmetric matrix, with components Cij 

being the smallest constant satisfying 

{TiV,Tjw)A < eij{TiV,v)\{TjW,w)X, \/v,w ^V. (3.14) 


Then we have 

Ki < p{e), 

where p{e) denotes the spectral radius of matrix e. Furthermore, if Cij < for some 7 G (0,1), 

then p{e) < (1 - 7 )"^ 

4 Convergence Analysis 

We here first introduce two interpolation norms and relevant Sobolev spaces (see e.g., [10]). Let A 
be a domain in R^. For integer m, denote by || • the Sobolev norm of integer order m, i.e., 


1/2 


I A’" (A) 




\<m 


with I = {li,l 2 ), \l\ = h + h and Let p > 0 he a non-integer and 0 < s < 1 , n is 

a non-negative integer such that n < p < n + 1. We introduce the interpolation norms 




f r°° _ \ / r°° ^ \ 

Vo ^W'^WnpA) ■= K{v,t)t~‘^‘'~^dtj (4.1) 


where 


K{v,t):= mf (||w+t ||i«||^„+i 


2||„..l|2 

H’^+i(A) J ’ 


K{v,t):= inf I|T-W||^ 2 (a) 

weH^(A) ^ 




Relevant Sobolev spaces are 

#^(A) := {u G l2(A);11u| 


(A) 


< 00 }, H^{A) := {v G L‘^(A)]\\v\ 


HpA) 


< 00 }. 


(4.2) 
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Let Ai, A 2 be two domains in with Ai C A 2 , and then 


/ inf 

/o 






< 


/ inf 

^Jo «;eA"+i(A2) 

7“ j.f (1 

So we have, for v G H^{A 2 ), 


^ ^)UillHn(Ai) ^ II''^|AilljJn+i(Ai) 1 ^ 


< 


V '“^IIj5'"(A2) * II'“^IIA"+1(A2) 


I^IIam(Ai) - II^IIam(A2)' 




1/2 


\ 1/2 

-‘^^^-^dt j 

1/2 


(4.3) 

(4.4) 


Remark 4.1. The following space relations can be found in literature: (1) /i > 0, iLi‘(R^) and Hq{Q) 
coincide with and LfQ(n) respectively; (2) for 1/2 < < 1, Hq{Q) coincides with H^{Q) 

(see [ 22 II 3 O]); for 1/2 < /t < 1, Hq{Q) coincides with Hq{Q,) (this can been shown by (1), (2) and 
the definitions of the interpolation spaces). 

Combining with remark 14.11 and the well known interpolation property (see e.g., Lemma 22.3 
in m), we know, for 1/2 < /i < 1, 


||(/-Qfc)ii|| </ifc||ii||//M(o), veH^{n). (4.5) 

Now, we develop some results for the finite element spaces 14, A: > 1. Let fl' C be a suitable 
polygonal domain such that Q C 0,' and dist((90', 14) > C for a positive C. Tj^, k > 1, are the 
quasi-uniform triangulations obtained by extending Tk from 14 to 14', that is, 7/' in 14 coincides 
with Tk- Furthermore we still make sure that 7/' = {t^} consists of simplexes of size hk- Let 
V/ = {u G C'^(i4') : = 0, u|t- G G 7/'}. In the following, for v G 14, v always denotes 

its extensions (on 14' and on M^), which is extended by zero outside 14, and so we also have u G V/. 

Lemma 4.2. Let > 0, v £ with supp{v) C 14 /u also denotes its extension on which 

is extended by zero outside 14'/. Then we have | bl |//m(o') ~ ll'f^llH''(R2)- 


Proof. For fj, being a integer, the conclusion is direct. For the case that pi is not a integer, denote n 
as a non-negative integer such that n < p, < n+1. From (14.311 . ||^^||//M(r 2 ') ^ I |'^’||j 7 m(r 2 ) ~ I |'^’||hm(r 2 ). 
Now we prove the converse relation. Let A be a domain in with smooth boundary such 

that 14 CC A C LI'. Then by (14.4p . v G H^{K). Following the proof for the strong extension 
of Sobolev space (see e.g.. Theorem 4.26 in m), we can show that there is a linear operator E 
continuous from H^{A) into iL^M^) for integers 0 < / < n -|- 1, such that £'(^4) = v. Then we 
have 
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00 

inf 
0 weH^+^ 


V — w 


w 


+ i7kll|n + l 




1/2 


< 


< 

rs_/ 


r 

y7o w&H^+^{A) 

(r 

y7o weH^+^{A) 

I^lli5-(*{A)> 


{\\E{v\A-w)\\%„^^,.^+t^\\Ew\\% 


2 _ 

JJn + l 


\ 1/2 

t~^f^~^dt I 




V '“'llHn(A) + A llll'IljJn+qA) 


^ t~^>^~^dt 


1/2 


(4.6) 
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where the last inequality is by the continuity of E. Combining with (14.41) . we obtain ||u||j:^m(]r 2 ) ~ 
Lemma 4.3. For 0 < fi < 3/2, v G 14, we have 

^ (4.7) 

and then 14 C 77^ (M^). 

Proof. By u G 14, we known v G Vf, from [3l li71l48] and further ||u||j;^m(r 2 ) < 

/i/'^||u|| bv Lemma 14.21 □ 

Let /3 be a positive with a + ft < 3/2 and a — /3 > 0 in the rest of this paper. We have the 
following results: 

Lemma 4.4. It holds that 


{v,w)a < ||^^||hc+/3(r2)||u;||j:^c-,3(r2), v,w G V. 
Proof. Since v,w ^V, by Lemma 14.31 we know that v,w ^ 77“^^(M^). Then 
{v,w)a = {Av,w) = B{v,w) 


p27V 

/ M{e)de+ c{v,w) 

Jo 

f2TT 

(77“+^, M{6)d6 + c{v, w) 


r27T 

/ \\Dg^^v\\L2(^^)\\Dg~^w\\L2(^Q)M{e)de + c || u |42( q ) || u ;|42( q ) 

Jo 

||u||^a+/3(R2)||u;||j^Q,-/3(]K2) + c||u|42(q) | |'fi'| | L2 (H) 
\\v\\j^a+/3(^^2\\\w\\fja-li 


< 

< 

rs_/ 

< 


(4.8) 


where the third equality is by Lemma 12.111 and the second inequality is by Lemma 12.101 and 

j^^M{e)de<i. □ 

Lemma 4.5. Let i < j, then 

{v,w)a < ||w^||, V e Vi, w e Vj. 

Here we recall that 7 G (0,1) is a constant such that h^ = 0{'y^). 

Proof. For v G Vi, w G Vj, we know that 

{v,w)a < ||^^||//c+/3(r2)||'«;||j;^c.-/3(r 2) < /i7("“^)||u;| 

= {hj/h,fhr»h-‘^\\v\\ ||u;|| < 


u re 


where the first inequality is by Lemma 14.41 the second inequality is by Lemma 14.31 and the last 
inequality is by the relation h^ ~ 0 ( 7 ^). □ 

Lemma 4.6. Let Wi = {Qi — Qi-i)V, then 


{v,w)A<'y^^ 4/3||^|m|y;m^ \Iu GWi,V GWj. 


(4.9) 
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Proof. By (14.511 and (|3.3I) . we have 


Ibll ^ ^fcll^lU, 'ivGWk- 

Combining the above inequality with Lemma 14.51 gives the lemma. 

Lemma 4.7. It holds that 

{TiV, Tjw)A < v)\{TjW, w)\, Vv, w G V. (4.10) 

Proof. It suffices to prove (j4.inp holds for i < j. Assume that i < j, and then for v,w G V, 


{TiV,Tjw)A = {RiAiPiV,RjAjPjw)A 

< jU-i)hh-<^h-°\\R^AiPiv\\\\RjA,Pjw\\, (4.11) 

where the inequality is by Lemma 14.51 

\\RiAiPiv\f = {RiAiPiV, RiAiPiv) « ^{RiAiPiV, AiPiv) = ^{TiV,v)A, 

where the second equality is by (13.7|) and the symmetry of R^. Then we obtain 

\\RiARv\\<Xi^^\TiV,v)][\ (4.12) 

and similarly 

\\RjAjP,w\\ < X-^^\TjW,w)f. (4.13) 

For u G 14, we have 

{Av,v) = \\v\\\ « ||w|Ih-(C!) ^ (4-14) 

where the second equality is by (13.3p and the last inequality is by Lemma 14.31 For w G V, 
V = (Qk - Qk-i)'w G 14, by (Ii3]) . we have 

~ {Av,v). (4.15) 

By (|4.141) and (14.15p . it is not hard to see that 

Xk-hf^, k = l,2,...,J. (4.16) 


Combining (14.1111 with (|4.12p . (|4.13p and ()4.16p gives 


{TiV,Tjw)A<7^^ "''^{TiV,v)l{TjW,w)X, '^v,wGV. 


The Lemma is proved. □ 
Lemma 4.8. Let 


and then for v gV , we have 


J 


'^Wm •— '^^ 11 {Qk Qk—l)v\\A'! 

k=l 


M\m 


TIU- 


(4.17) 
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Proof. It is not hard to see that the space coincides with in [30]. Combining with 

Theorem 1 of [30], we know that ~ “ Qk-i)w\\^ holds for w G 

For V e V, \\{Qk - Qk-i)v\\^ ~ - Qk-i)v\\\jc^Q) by (|4.14p and (I4.15p . Combining with 

(13.3p gives the lemma. □ 

Theorem 4.9. ITe have 

Ko<l, 

That is to say, our V-cycle multigrid method is optimal, which means that the convergence rate is 
independent of the mesh size and mesh level. 

Proof. For v G V, decompose v as v = Yli=i with Vi = {Qi — Qi-i)v. By (|4.5I) and (|3.3I) we have 
lk*ll ^ Furthermore combining (13.7p with (14.161) . we have {R~^Vi,Vi) < ||uj||^. Using 

Lemma 14.81 gives Kq < 1. Finally Combining Lemma 14.71 with Lemma 13.41 gives that Ki <1. □ 


5 Implementation 

For simplicity, in this section, we only consider / = 1 in p3.4p . i.e., 14, k = 1,2,..., J, are the spaces 
consisting of the piecewise linear polynomials. Let f)],, i = 1,..., be the nodal basis of the 
finite element space 14. The implementation are a classical procedure in literature (see e.g., 121 ), 
and we here only illustrate how to generate the stiff matrices of the finite element systems and how 
to choose Rk '. Vk ^ Vk, k = 2,..., J, the approximations of A^. 


5.1 The stiffness matrices and 

For Ak, denote its corresponding stiffness matrix by A^ G with entries 

(5.1) 

Since M has the discrete form M{9) = ~ ^0, 


{^k)ij 


L 


Y, Pi{DlYf^k.De,+.fP^) + c( 4 , -^i). 

1=1 


We need only discuss how to numerically compute 

= [ Dg'^Dgcfl X Dg+^cfldxdy (5.2) 

Jssup{ 4 >i) 

for a hxed 6, where n = (2 — 2a), and then the entries of the stiff matrices can be numerically 
computed. If a = 1, the computation of the stiffness matrices is easy, since the original problem is 
an integer order one. Now we focus on the case of 1/2 < a < 1. Define the index set Ki as 

Ki = {l-,Ti GTk,Tic supp(0y}. 


Then 

R = Y [ ^e^'^d^k X Dg+.,,(l){dxdy 

leKj 

l£Kj I'GKi 
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Figure 1: Illustration for computing Iq 


where for a set S in 


Xs{x,y) = 


1, if {x,y) G S’; 


0 , otherwise. 

Noting that D 0 {<fl.)\^ii, are both constants, we numerically compute 


4 


Dg '^x.^v (x, y) X y I (x, y)dxdy, 


(5.3) 


and then Ig can be computed. Next we illustrate how to compute the integral in (|5.3p by an example. 
On the left of Figured] is Cartesian coordinate systems xOy and x'Oy', and the angle between axes 
Ox and Ox' is 9. On the right of Figured! the two triangles are and r^; Di, D 2 ,D^ denote the 
corresponding vertices of the triangles; fl/, 0 ,jj denote the corresponding shadow areas respectively; 
lines DiPi and D 3 P 2 are both Parallel to axis Oy'] 71,72,73,74,75,76 are correspondent angles. 
Denote the coordinates of Di,D 2 and D 3 under coordinate system x'Oy' by {x'ii,y'i 2 )-, (a^ 2 i) 2 / 22 ) 
and ( 3 ^ 31 , 2 / 32 ) respectively. Then we have 


Dg {x, y) X X / (x, y)dxdy 


/ Dg''Xri'ix,y)dxdy + / Dg''x.^i'ix,y)dxdy 

iQi * JQii 


- 00 -D {x', y')dx'dy' + / _ooD Tx ,/ (x', y')dx'dy' 




1 


(x' - x'l + {y'l - y') ianxiY dx'dy' 


Inr r(z^ + 1) 


1 


(x' - x'l - (y'l - y') tan72)'' dx'dy' 


ni r(z^ +1) 

(x' - x'2 + {y' - ya) tan 73 )"" dx'dy' 


Inn + 1 ) 


1 


r(i. + l) 


(x' - X 2 + {y' - y'2) tan 74 )'' dx'dy'. 


The last four integrals above can be computed directly. Finally we know that the entries of the 
stiffness matrices can be numerically computed. 

We choose as 

Nk 


RkV = V € 14, 

Afc 


(5.4) 
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Figure 2: Uniform triangulation 


with Afc ~ ^kh\- Define mass matrix E with entries 

{Mkh = 

For V E Vfc, denote v E the vector of coefficients of v in the basis It is known that 

v'^MkV ~ h^v'^v and %P"M'^v ~ h\v^M^v. Hence we have 

{RkV, v) « ^v^mIv « ^v^MkV ^ -^(u, v), (5.5) 

Afc ‘^k ■^k 

which means (13.7p holds. In the numerical tests, we take A^ = A: = 2,..., J. It is not hard 

to verify that {Ak)ii ^ h\\k. 


5.2 Computation complexity 

For the numerical approximation of SFPDEs, one of the key issues is how to reduce the computation 
complexity. We confine ourself to the case that D is a square domain, and of course the technique 
here is also helpful for effectively designing schemes for the case that D is a general domain. 

The triangulations 7fc, k = 1,2,...,J are those in Figure [21 where dashed curve denote the 
ellipsis, Uk = no2^ — 1, /^ = lo2^ — 1 with positive integers no, lo, and p™, m = 1,..., Ukh are the 
interior points. The finite element space 14 = {u E Hq{Q,) : v\r E Pi(t),Vt E Tk}- Let 0™ = 
m = 1,... ,nklk, be the nodal basis functions, i.e., 0™ is a piecewise linear polynomial 
whose values are 1 at p™ and zeros at other nodes (including interior and exterior nodes). 

Denote U = (Ui, U 2 ■ ■ ■, Un^., ■ ■ ■, U2nk,: • • • j Next we discuss how to effectively con¬ 

duct the multiplication of matrix Ak and vector U. Let u = {i'i,V2, ■ ■ ■ ,i^(2nk-i)ik-nk+i)'^ ^ 

—l)/fe—^fc+l 


n( 2 nfe-i)j+i = H((/.fc, i = 1,... ,nfc,j = 0,... ,4 - 1, 

iy( 2 nk-i)j-i +2 = z = 2,... ,nfc, j = 1,... ,4 - 1- 

Define a symmetric Toeplitz matrix 


/ 


A = 


121 

1^2 


^2 • • • l'(^2nk-l}lk-nk ^(2nk-l)lk-nk+l ^ 

^1 * * * ^(2nfc —!)//(.—nfc — 1 ^(2?!^ —1)/^— 


\ ^{2nk-l)lk-nk+l ^(2nk-l)lk-nk 
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122 


122 

121 


/ 















































Toeplitz matrix, also called diagonal-constant matrix, is a matrix in which each descending diagonal 
from left to right is a constant. 

For any i, j with 1 < i < j < > let dj, r*, dj , rj be nonnegative integers satisfying i = Ukdi+ri, 

j = Ukdj + Tj, 1 < ri,rj < Uk- Let j' = {dj — di),i' = \rj — ri\, and then by the property of the 
operator B{-,-), it is easy to see that 


B{<i>l4) 


if rj > n; 
= z^/{ 2 nfc-i)-i'+i, if rj < n. 


And thereby any component of matrix A is also one of vector v. Define sets 


Im = {m{2nk - 1) -M, m{2nk - 1) -b 2,, m{2nk - 1) + rik}, m = 0,1,..., 4 - 1 

and I = U Im- We have the relation 


Afc = Ax^x, (5.6) 

where Ax^x denotes the sub-matrix of A which consists of entries Aij of A indexed by i,j E I. 
Denote U' E RC^^k-m-nk+i 

rife —1 

U (t/i 5 ■ ■ ■ 5 , 0, . . . , 0, ; • • • ? ; 0, . . . , 0, J72nfeH-l ? • ■ ■ ; ) * 

It is not hard to see that 

AkU = {Au')x, 

where for a given vector v, vx denotes the vector which consists of entries Vi indexed by z E X. So 
the multiplication of the matrix Ak and any vector U E can be obtained by conducting the 

multiplication of the Toeplitz matrix A and U' E R^‘^'^k-^)ik-nk+^_ multiplication of a Toeplitz 
matrix in and a vector in M”' can be done with computation complexity 0(nlogn). Recall 

that Nj = njlj denotes the number of the unknowns in the finite element problem (13.211 . and then 
by the above analysis, we conclude that for the V-cycle multigrid methods developed in Section 01 
each iteration needs computation complexity 0{Nj log Nj). 


5.3 Numerical results 

In this section, we shall present some numerical results to confirm our theoretical findings. In our 
numerical test, we take uq = Iq = 4, and take Nq = 4(nj -|- 1) if M is a continuous function. 

We shall check our V-cycle multigrid method and the preconditioned conjugate gradient algo¬ 
rithm (PCG) with Bj as the preconditioner. Meanwhile, the numerical result for the conjugate 
gradient algorithm (CG) is also presented for comparison. Our tests are carried out using Matlab 
software. The stopping criterion of the algorithm is 

\\u^ - u’^~^\\oo < 10 -®. 

We present two examples: one is with the probability measure M having a discrete form and 
the other with M being a continuous function. Table [T] and Table [2] list the numerical results 
for Example 15.11 and Example 15.21 respectively, where ”DOEs” denotes the degree of freedoms and 
’’Iter” denotes the iterative steps on each level. It is seen that the numbers of iterations of our 
V-cycle multigrid and PCG per level are bounded independent of the mesh size and mesh level, 
which confirms our theoretical results. 
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Figure 3: The CPU time per iteration 


Example 5.1. Let U = [0, 2] x [0, 2], the equation to be solved is 


1 

4* 


_ 1( ni-5 f n 

, I—OO-^-r ' \ —OO^ 


1.5 




(5.7) 


Level 

J 

DOFs 

V-cycle 

PCG 

CG 

Iter 

Iter 

Iter 

4 

4096 

13 

7 

58 

5 

16384 

13 

6 

72 

6 

65536 

13 

7 

118 

7 

262144 

13 

7 

197 

8 

1048576 

12 

7 

313 


Table 1: Numerical results for Example 15.11 


Example 5.2. Let U = [0,2] x [0,2] and M{6) = 1. The equation to be solved is 


- = 1. (5.8) 

We choose smooth f{x,y) in the examples such that the solutions have singularity near the 
boundaries. The computation complexity of our multigrid methods are shown in figure [3l where 
’’Time” denotes the CPU time (in seconds) spent by one iteration. As can be seen from the figure 
[31 the CPU time of each iteration ia almost linear with respect to the degree of freedoms. So the 
computation complexity of our multigrid method is also optimal. 


Level 

DOFs 

V-cycle 

PGG 

GG 

Iter 

Iter 

Iter 

4 

4096 

11 

6 

50 

5 

16384 

11 

6 

69 

6 

65536 

11 

6 

115 

7 

262144 

11 

6 

184 

8 

1048576 

11 

6 

307 


Table 2: Numerical results for Example 15.21 
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Appendix 

The Fourier analysis plays critical roles in this paper: for g G the Fourier transform of g is 

the function Tg defined on (the dual of) by 

-^5(6,6)= [ e~^^^^^^^^y^^'>g{x,y)dxdy, 

JR2 

where i denotes the imaginary unit; for g G the Fourier transform Tg of g is defined in the 

following distribution sense (see |4Uj 1: 

{Fg,v) = {g,Fv), Vu G Co“(m2), 

and more precisely, F is an isometry from L^(M^) into itself, which satisfies Parseval’s formula 
(see [21]) 

11-^511 = lbll 

and 

(u, uJ) = {Fv, Fw), 

where T denotes the complex conjugate of the complex number 2 . The Fourier transform of the /rth 
order fractional derivative consists of the complex in the form (zk)^ with /z > 0, k G M (see [32] 1. So 
it may be a multi-valued function. To guarantee the Fourier transform to be univalent, we express 
complex variable 2 ; = | 2 ;| exp(z0), —tt < 6 < tt, where exp(z0) = cos 9 + i sin 9, \z\ and 9 respectively 
denote the modulus and the argument of z. Then 

(zk)^ = (sign(K)z|/t|)^ = (|k| exp(zsign(K)7r/2))^ = |k|^ exp(z/zsign(K)7r/2), 

{— ik )^ = |k|^ exp(—z/zsign(K)7r/2). 

It is easy to see that, for ^ > 0, 


{—inY = {if^Y, V/t G M. 


(AT) 
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